########################################################
## Replication of SI
########################################################
rm(list=ls())
start_time <- Sys.time()
########################################################
## Package loading
########################################################
library(BridgeChange)
library(devtools)  
library(coda)
library(plm)
library(pforeach)
library(dplyr)
library(colorspace)
library(monomvn)
library(glmnet)
library(MASS)
library(MCMCpack)
library(pcse)
library(pastecs)
library(xtable)
library(bayesplot)
library(dotwhisker)
library(broom)
library(stableGR)
library(tidyverse)
library(lmtest)  

## set working directory for downloaded folder
## setwd("~/Dropbox/PAreplication/")
source("helper.R")
## setwd("~/GitHub/bayesianbridge/codes/ReplicationPA/")
dir.create("SI")

########################################################
## Figure 1 and Figure 2
########################################################

## We construct 24 sets of high dimensional TSCS data with varied group sizes ($n = (5, 10)$), time lengths ($t = (20, 40)$),  predictor sizes ($k = (50, 100)$), and break numbers ($m = (0, 1, 2)$) to evaluate the validity of our proposed method under a more challenging setting. 

########################################################
## Simulation Setting
########################################################
set.seed(39212)
n_chain <- 3
mcmc <- 1000
burn <- 1000
verbose <- 500
thin <- 1

t_time   <- c(20, 40)
n_unit   <- c(5, 10)
k_cov    <- c(50, 100)
n_break  <- c(0, 1, 2)
s_sparse <- c(1)

## 24 variations
sim_set <- expand.grid(n_unit, t_time, k_cov, n_break, s_sparse)
colnames(sim_set) <- c("N","T","K","B","S")

## sim_set <- expand.grid(n_unit, t_time, k_cov, n_break, s_sparse)
n.sim <- nrow(sim_set)

bias.report <- rmse.report <- matrix(NA, n.sim, 4) ## c(bias.ols, bias0, bias1, bias2)
gr.report <- waic.report <- matrix(NA, n.sim, 3) ## c(bias.ols, bias0, bias1, bias2)
state.report <-  matrix(NA, n.sim, max(t_time))

fixed = TRUE
########################################################
##  Simulation Starts!
########################################################
sim.name = paste0(lubridate::today(), "_sqrt_smallpanel_k", k_cov[1], "_n", n_unit[1])
dir.create(file.path(paste0("SI")))

for(cl in 1:nrow(sim_set)){

    N <- sim_set[cl, 1] 
    T <- sim_set[cl, 2] 
    K <- sim_set[cl, 3] 
    B <- sim_set[cl, 4]  
    S <- sim_set[cl, 5]
    Q <- 2; ## for random effects

    cat("\n==================================\n")
    cat("simulation = ", cl, "\n")
    cat("setup (T, N, K, B, S) = ", unlist(sim_set[cl,]), "\n")
    cat("==================================\n")
    
    ## generate data ==================
    NT <- N * T; 
    break.point <- round(T / 2)
    if(B == 2){
        break.point <- c(round(T/3), T - round(T/3))
    }
    break.sigma <- 1; 
    
    ## generate covariates
    x <- as.matrix(scale(abs(gen_cov(N, T, K)))) ## NT by K+1
    K_full <- ncol(x)
    
    ## ground truth
    mu1 <- 2
    K.ground <- ifelse(fixed, K, K + 1) ## if random, intercept is added.
    if (B == 0) {
        true.beta <- gen_beta(K.ground, mu = mu1, pattern = S, n_break = B)
        true.sigma <- sqrt(2)
        true.D     <- diag(runif(Q), Q)
        
    }else if (B == 1) {
        out <- gen_beta(K.ground, mu = mu1, pattern = S, n_break = B)
        
        true.beta1   <-  out[[1]]
        true.beta2   <-  out[[2]]
        true.beta <- list(true.beta1, true.beta2)
        true.sigma <- list(sqrt(2), sqrt(2));
        true.D     <- list(diag(runif(Q), Q), diag(runif(Q), Q));
 
    } else {
        out <- gen_beta(K.ground, mu = mu1, pattern = S, n_break = B)
        true.beta1   <-  out[[1]]
        true.beta2   <-  out[[2]]
        true.beta3   <-  out[[3]]
        true.beta <- list(true.beta1, true.beta2, true.beta3)
        true.sigma <- list(sqrt(2), sqrt(2), sqrt(2));
        true.D     <- list(diag(runif(Q), Q), diag(runif(Q), Q), diag(runif(Q), Q));
        
    }
    
    epsilon    <- 0
    
    
    
########################################################
    ## DGP
########################################################
    ## x <- gen_cov(N, T, K)
    out <- dgp.break(x, N = N, T = T, Q = Q, B = B,
                     true.beta, true.sigma, true.D,
                     break.point = break.point,  break.sigma = 1,
                     fixed=TRUE,  print = FALSE )

    ## true state and beta
    
    
    ## raw data
    y <- out$y
    X <- out$X
    subject.id <- out$subject.id;
    time.id <- out$time.id
    time.dummy <- out$time.dummy
    data = as.data.frame(cbind(y, X, subject.id, time.id))
    true.beta.time <- out$true.beta.time
    coplot(y ~ time.id|subject.id, data=data)

    
########################################################
    ##  Estimation
########################################################

    ## effect = "twoways": do not use it
    effect = "individual"
    model = "within"
    W <- matrix(0, length(y), 1)
    
    measurevar <- "y"
    groupvars  <- colnames(X)
    
    ## This creates the appropriate string:
    formula <- as.formula(paste(measurevar, paste(groupvars, collapse=" + "), sep=" ~ "))
    index = c('subject.id', 'time.id')
    ## model.matrix <- getFromNamespace("model.matrix.pFormula", "plm")
    pdata    <- plm::pdata.frame(data, index)
    ## pformula <- plm::pFormula(formula)
    
    plm.X <- model.matrix(formula, pdata, rhs = 1, model = model, effect = effect)
    plm.y <- pmodel.response(formula, pdata, model = model, effect = effect)
    effect = "individual"
    model = "within"
    ols <- plm(formula, pdata, model = model, effect = effect)
    ols.beta <- ols$coef
    
    unscaled.Y <- plm.y
    unscaled.X <- plm.X
    subject.id <- as.numeric(as.factor(data[,index[1]]))
    time.id    <- as.numeric(as.factor(data[,index[2]]))
    
    c0 = .1; d0 = .1; r0 =  1; R0 = 1;
    
    sim_out <- pforeach(i = 1:n_chain, .cores = n_chain, .seed = 3821, .c = "list")({
        test0 <- BridgeMixedPanel(subject.id = subject.id, time.id = time.id,
                                  standardize = FALSE, ## beta.alg = "BCK", 
                                  mcmc = mcmc, burn = burn, thin = thin, verbose = verbose,
                                  y = y, X = X, W = W, n.break = 0, 
                                  r0 = r0, R0 = R0, c0=c0, d0=d0,
                                  alpha.MH = TRUE,
                                  unscaled.Y= y, unscaled.X = X, Waic = TRUE)
         
        test1 <- BridgeMixedPanel(subject.id = subject.id, time.id = time.id,
                                  standardize = FALSE, 
                                  mcmc = mcmc, burn = burn, thin = thin, verbose = verbose,
                                  y = y, X = X, W = W, n.break = 1, 
                                  r0 = r0, R0 = R0, c0=c0, d0=d0,
                                  alpha.MH = TRUE,
                                  unscaled.Y= y, unscaled.X = X, Waic = TRUE)
         
        test2 <- BridgeMixedPanel(subject.id = subject.id, time.id = time.id,
                                  standardize = FALSE, 
                                  mcmc = mcmc, burn = burn, thin = thin, verbose = verbose,
                                  y = y, X = X, W = W, n.break = 2, 
                                  r0 = r0, R0 = R0, c0=c0, d0=d0,
                                  alpha.MH = TRUE,
                                  unscaled.Y= y, unscaled.X = X, Waic = TRUE)
        list("test0" = test0, "test1" = test1, "test2" = test2)
    })  
    
    ## stack the result
    test0l <- mcmc.list(lapply(sim_out, function(x) x[[1]][,grep("beta", colnames(x[[1]]))]))
    test1l <- mcmc.list(lapply(sim_out, function(x) x[[2]][,grep("beta", colnames(x[[2]]))]))
    test2l <- mcmc.list(lapply(sim_out, function(x) x[[3]][,grep("beta", colnames(x[[3]]))]))
    
    ## further stack the result (from multiple chains)
    test0 <- do.call("rbind", test0l)
    test1 <- do.call("rbind", test1l)
    test2 <- do.call("rbind", test2l)
    
    est.beta0 <- test0[, grep("beta", colnames(test0))]
    est.beta1 <- test1[, grep("beta", colnames(test1))]
    est.beta2 <- test2[, grep("beta", colnames(test2))]
    
    ## m = 0
    mat.beta0 <- matrix(rep(apply(est.beta0, 2, mean), T), dim(est.beta0)[2], T)
    
    ## m = 1
    mu.beta1 <- apply(est.beta1, 2, mean)
    mu.beta11 <- mu.beta1[grep("regime1", names(mu.beta1))]
    mu.beta12 <- mu.beta1[grep("regime2", names(mu.beta1))]
    state1 <- round(apply(attr(sim_out[[3]][[2]], "s.store"), 2, mean))
    mat.beta1 <- matrix(rep(mu.beta11, T), dim(est.beta0)[2], T)
    mat.beta1[, state1==2] <- mu.beta12
    
    ## m = 2
    mu.beta2 <- apply(est.beta2, 2, mean)
    mu.beta21 <- mu.beta2[grep("regime1", names(mu.beta2))]
    mu.beta22 <- mu.beta2[grep("regime2", names(mu.beta2))]
    mu.beta23 <- mu.beta2[grep("regime3", names(mu.beta2))]
    state2 <- round(apply(attr(sim_out[[3]][[3]], "s.store"), 2, mean))
    mat.beta2 <- matrix(rep(mu.beta11, T), dim(est.beta0)[2], T)
    mat.beta2[, state2==2] <- mu.beta22
    mat.beta2[, state2==3] <- mu.beta23
    
    ## plm estimates: twoway fe
    ## effect = "twoways": do not use it
    mat.ols <- matrix(rep(ols.beta, T), dim(est.beta0)[2], T)
    
    
    bias0 <- mean(true.beta.time - mat.beta0)
    bias1 <- mean(true.beta.time - mat.beta1)
    bias2 <- mean(true.beta.time - mat.beta2)
    bias.ols <- mean(true.beta.time - mat.ols)
    
    rmse0 <- rmse.fn(true.beta.time - mat.beta0)
    rmse1 <- rmse.fn(true.beta.time - mat.beta1)
    rmse2 <- rmse.fn(true.beta.time - mat.beta2)
    rmse.ols <- rmse.fn(true.beta.time - mat.ols)
    
    ## Approximate convergence is diagnosed when the upper limit is close      to 1
    ## Potential scale reduction factors:
    gr <- c(stableGR::stable.GR(test0l)$mpsrf, stableGR::stable.GR(test1l)$mpsrf, stableGR::stable.GR(test2l)$mpsrf)
    
    ## state.report
    rmse.report[cl, ] <- c(rmse.ols, rmse0, rmse1, rmse2)
    gr.report[cl, ] <- gr
    if(B==0){
        
    }else if(B == 1){
        state.report[cl, 1:length(state1)] <- state1
    }else{            
        state.report[cl, 1:length(state2)] <- state2
    }
    
    ## waic list
    waic.report[cl, ]<- c(mean(unlist(lapply(sim_out, function(x) {attr(x[[1]],  "Waic.out")[1]}))),
                          mean(unlist(lapply(sim_out, function(x) {attr(x[[2]],  "Waic.out")[1]}))),
                          mean(unlist(lapply(sim_out, function(x) {attr(x[[3]],  "Waic.out")[1]}))))
    
    
    ## save file here
    saveRDS(rmse.report, file = paste0("SI/rmse.rds"))
    saveRDS(gr.report, file = paste0("SI/gr.rds"))
    saveRDS(waic.report, file = paste0("SI/waic.rds"))
    saveRDS(state.report, file = paste0("SI/state.rds"))
    
########################################################
    ## Step 6. Store
########################################################
    
    cat("\n==================================\n")
    cat("simulation = ", cl, " is done ! \n")
    cat("==================================\n")
    
}

rmse <- readRDS(file = paste0("SI/rmse.rds"))
gr <- readRDS(file = paste0("SI/gr.rds"))
waic <- readRDS(file = paste0("SI/waic.rds"))
state <- readRDS(file = paste0("SI/state.rds"))
 

##
true.m <- sim_set$B
true.t <- sim_set$T
true.n <- sim_set$N
true.k <- sim_set$K
true.state <- list()
for(i in 1:n.sim){
    if(true.m[i] == 1){
        true.state[[i]] <- c(rep(1,  true.t[i]/2), rep(2,  true.t[i]/2))
    }else if(true.m[i] == 2){
        break.point <- c(round(true.t[i]/3), true.t[i]- round(true.t[i]/3))
        true.state[[i]] <- c(rep(1,  break.point[1]), rep(2,  break.point[2] - break.point[1]),
                             rep(3,  true.t[i] - break.point[2]))
    }else{
     cat("m = ", true.m[i], "\n")   
    }
}


## Graph 1: hidden state
pdf(file = paste0("SI/Figure2a.pdf"),
    width=12, height = 8, family="sans")
par (mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.01)
par(mfrow=c(2, 2))
## subplot 1
plot(1:max(true.t[true.t==unique(true.t)[1] & true.m == 1]), ylim=c(0.5, 2.5), type="n", axes=FALSE, ylab="Hidden state", xlab="Time");
axis(1); axis(2, at=1:3); box();mtext(paste0("Hidden State Recovery (m = 1, T = ", unique(true.t)[1], ")"), 3)
for(i in which(true.t==unique(true.t)[1]& true.m == 1)){
    lines(jitter(state[i, 1:length(true.state[[i]])]), col="grey60", cex=0.2)
    lines(true.state[[i]], col="grey20", lwd=2)}
legend("topleft", legend=c("True", "Estimated"), col=c("grey20", "grey60"),
       lty=1, lwd=c(2, 1), bty="n")

## subplot 2
plot(1:max(true.t[true.t==unique(true.t)[1] & true.m == 2]), ylim=c(0.5, 3.5), type="n", axes=FALSE, ylab="Hidden state", xlab="Time");
axis(1); axis(2, at=1:3); box();mtext(paste0("Hidden State Recovery (m = 2, T = ", unique(true.t)[1], ")"), 3)
for(i in which(true.t==unique(true.t)[1]& true.m == 2)){
    lines(jitter(state[i, 1:length(true.state[[i]])]), col="grey60", cex=0.2)
    lines(true.state[[i]], col="grey20", lwd=2)}
legend("topleft", legend=c("True", "Estimated"), col=c("grey20", "grey60"),
       lty=1, lwd=c(2, 1), bty="n")

## subplot 3
plot(1:max(true.t[true.t==unique(true.t)[2] & true.m == 1]), ylim=c(0.5, 2.5), type="n", axes=FALSE, ylab="Hidden state", xlab="Time");
axis(1); axis(2, at=1:3); box();mtext(paste0("Hidden State Recovery (m = 1, T = ", unique(true.t)[2], ")"), 3)

for(i in which(true.t==unique(true.t)[2]& true.m == 1)){
    lines(jitter(state[i, 1:length(true.state[[i]])]), col="grey60", cex=0.2)
    lines(true.state[[i]], col="grey20", lwd=2)}
legend("topleft", legend=c("True", "Estimated"), col=c("grey20", "grey60"),
       lty=1, lwd=c(2, 1), bty="n")

## subplot 4
plot(1:max(true.t[true.t==unique(true.t)[2] & true.m == 2]), ylim=c(0.5, 3.5),  type="n", axes=FALSE, ylab="Hidden state", xlab="Time");
axis(1); axis(2, at=1:3); box();mtext(paste0("Hidden State Recovery (m = 2, T = ", unique(true.t)[2], ")"), 3)
for(i in which(true.t==unique(true.t)[2]& true.m == 2)){
    lines(jitter(state[i, 1:length(true.state[[i]])]), col="grey60", cex=0.2)
    lines(true.state[[i]], col="grey20", lwd=2)}
legend("topleft", legend=c("True", "Estimated"), col=c("grey20", "grey60"),
       lty=1, lwd=c(2, 1), bty="n")
dev.off()


## Graph 2: diagnostics
pdf(file = paste0("SI/Figure2b.pdf"),
    width=8, height = 6, family="sans")
par (mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.01)
plot(sapply(1:nrow(sim_set), function(i)gr[i, true.m[i]+1]), pch=19,
     type="p", axes=FALSE, ylab="Rhat", xlab="Simulation")
axis(1); axis(2); abline(h=c(1.0, 1.1), lty=3, col="grey60");
box();mtext("Gelman-Rubin Convergence Diagnostics", 3)
dev.off()


## Graph 3: waic
pdf(file = paste0("SI/Figure1b.pdf"),
    width=15, height = 5, family="sans")
par (mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.01)
par(mfrow=c(1, 3))
counter <- 1
for(i in 1:3){##
    df <- waic[true.m == i-1, ]
    plot(df[1,], ylim = range(df), pch=19, axes=FALSE, type="n", ylab="WAIC", xlab="Model")
    for(j in 1:nrow(df)){
        lines(df[j,], col="grey70")
        points(df[j,], pch=19, cex=1)
        points(true.m[counter]+1, df[j, true.m[counter]+1], pch=19, col=NetworkChange:::addTrans("brown", 100), cex=5)
        text(true.m[counter]+1, df[j, true.m[counter]+1], paste0("T=", true.t[counter], ", N=", true.n[counter], ", K=", true.k[counter]),
             pos = 3 - (true.m[counter] - 1))
        counter <- counter + 1
        cat("counter  = ", counter, "\n")
    }
    mtext(paste0("Break number = ", true.m[counter-1]), 3)
    axis(1, at = 1:3, labels=c("break 0","break 1","break 2")); axis(2);  box()
}
dev.off()
    
    

    
## Graph 4: rmse
pdf(file = paste0("SI/Figure1a.pdf"),
    width=15, height = 5, family="sans")
par (mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.01); par(mfrow=c(1, 3))
counter <- 1
for(i in 1:3){##
    df <- rmse[true.m == i-1, ]
    plot(df[1, ], ylim = range(df), pch=19, axes=FALSE, type="n", ylab="RMSE", xlab="Model")
    for(j in 1:nrow(df)){
        lines(df[j,], col="grey70")
        points(df[j,], pch=19, cex=1)
        if(true.m[counter] == 0){
            ## plm
            points(true.m[counter]+1, df[j, true.m[counter]+1],
                   pch=19, col=NetworkChange:::addTrans("brown", 100), cex=5)
        }
        ## break 
        points(true.m[counter]+2, df[j, true.m[counter]+2],
               pch=19, col=NetworkChange:::addTrans("brown", 100), cex=5)
        text(true.m[counter]+2, df[j, true.m[counter]+2],
             paste0("T=", true.t[counter], ", N=", true.n[counter], ", K=", true.k[counter]),
             pos = 3 - (true.m[counter] - 1))
        counter <- counter + 1
        cat("counter  = ", counter, "\n")
    }
    mtext(paste0("Break number = ", true.m[counter-1]), 3)
    axis(1, at = 1:ncol(df), labels=c("FE", "HMBB break 0","HMBB break 1","HMBB break 2")); axis(2);  box()
}
dev.off()

########################################################
## Figure 3: HMBB_sim.pdf
########################################################
set.seed(11199)
## simulate data
K <- 80
n <- 100
X <- matrix(rnorm(n*K), n, K)
sig2 <- 4
beta.true <- matrix(NA, 2, K)

beta.true[1,] <- matrix(rnorm(K, 1, 1), K, 1)*2
beta.true[2,] <- matrix(rnorm(K, -1, 1), K, 1)
# beta.true[1, sample(1:K, K/2, replace=FALSE)] <- 0
# beta.true[2, sample(1:K, K/2, replace=FALSE)] <- 0
mu1 <- X[1:(n/2), ]%*%beta.true[1,]
mu2 <- X[((n/2)+1):n, ]%*%beta.true[2,] 
Y  <- c(rnorm(n/2, mu1, sqrt(sig2)), rnorm(n/2, mu2, sqrt(sig2)))

## fit he model
formula <- Y ~ X
fit.cp0 <- BridgeChangeReg(formula=formula, mcmc=1000, burn=1000, n.break = 0, Waic = TRUE)
fit.cp1 <- BridgeChangeReg(formula=formula, mcmc=1000, burn=1000, n.break = 1, Waic = TRUE)
fit.cp2 <- BridgeChangeReg(formula=formula, mcmc=1000, burn=1000, n.break = 2, Waic = TRUE)


## model selection by WAIC
## WAIC prefers the no break model but the ground truth is the one break model
(waic <- WaicCompare(list(fit.cp0, fit.cp1, fit.cp2), print = TRUE))
## obtain beta.hat
beta.est <-  matrix(apply(fit.cp1[, grep("beta", colnames(fit.cp1))], 2, mean), 2, , byrow=TRUE)

########################################################
## Figure 3
########################################################

## plot generated data
pdf(file = paste0("SI/Figure3.pdf"),
    width=12, height = 3, family="sans")
par (mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.01); 
par(mfrow=c(1,4))

## waic
plotWaic(waic)

## regime 1
plot(beta.true[1,], beta.est[1,],  xlab="TRUE", ylab="EST", type = "n",main="Regime 1",
       xlim = range(beta.est), ylim = range(beta.true), asp = 1)
abline(a=0, b=1, col="red", lty = 3, lwd = 1.5)
points(beta.true[1,], beta.est[1,], col="darkblue")

## regime 2
plot(beta.true[2,], beta.est[2,],  xlab="TRUE", ylab="EST", type = "n", main="Regime 2",
       xlim = range(beta.est), ylim = range(beta.true), asp = 1)
abline(a=0, b=1, col="red", lty = 3, lwd = 1.5)
points(beta.true[2,], beta.est[2,], col="darkblue")

## all covariatess
dotplotRegime(fit.cp1, hybrid=FALSE, location.bar=12, x.location="default",
              text.cex=0.8, main="Time-varying Movements of All Covariates")
dev.off()

########################################################
## Figure 4
########################################################

## true break at t = 50
pdf(file = paste0("SI/Figure4.pdf"),
    width=10, height = 5, family="sans")
par (mar=c(3,3,2,1), mgp=c(2,.7,0), tck=-.01); 
par(mfrow=c(1, 2))
MCMCpack::plotState(fit.cp1, legend.control =c(60, 0.85), main="One break")
MCMCpack::plotState(fit.cp2, legend.control =c(60, 0.85), main="Two breaks")
dev.off()



########################################################
## Big Simulation
########################################################
library(faux)
require(BridgeChange)
require(glmnet)
require(monomvn)
require(doParallel)
require(doMC)
require(doRNG)
require(mvnfast)
coef_bridge <- function(out){
    coef <- apply(out1, 2, mean)
    coef[grepl("beta", names(coef))]
}

## some parameters
pdf.options(family = "sans")
length.out <<- 5
delta <- rho <- seq(0.1, 1, length.out = length.out)
p      <- 200 # this is fixed
sigma  <- 2
n_iter <- 20
n_miss <- 0.2

## y: ceiling(n * rho): number of sparse signal
## 10  58 105 153 200
## x: floor(p * delta): number of observations
## 10  57 105 152 200

## Case 1:  No correlation no changepoint
source("nochange/corr00/pred_mse_sim_cluster_corr00.R")
source("nochange/corr00/coef_mse_sim_cluster_corr00.R")
source("nochange/corr00/cv_mse_sim_cluster_corr00.R")

##
## Table 3 No correlation no changepoint
##
source("nochange/corr00/pred_mse_sim_summary.R")
source("nochange/corr00/coef_mse_sim_summary.R")
source("nochange/corr00/cv_mse_sim_summary.R")

##
## Figure 5 No correlation no changepoint
##
source("nochange/corr00/pred_mse_sim_plot.R")
source("nochange/corr00/coef_mse_sim_plot.R")
source("nochange/corr00/cv_mse_sim_plot.R")


##
## Case 2: Correlated no changepoint simuuation
##
source("nochange/corr07/pred_mse_sim_cluster_corr07.R")
source("nochange/corr07/coef_mse_sim_cluster_corr07.R")
source("nochange/corr07/cv_mse_sim_cluster_corr07.R")

##
## Table 4 Correlated no changepoint
##
source("nochange/corr07/pred_mse_sim_summary.R")
source("nochange/corr07/coef_mse_sim_summary.R")
source("nochange/corr07/cv_mse_sim_summary.R")


##
## Figure 6 Correlation no changepoint
##

source("nochange/corr07/pred_mse_sim_plot_corr07.R")
source("nochange/corr07/coef_mse_sim_plot_corr07.R")
source("nochange/corr07/cv_mse_sim_plot_corr07.R")




##
## Case 3: correlation single changepoint
##
source("change/corr07/pred/pred_mse_change_sim_corr07.R")
source("change/corr07/coef/coef_mse_change_sim_corr07.R")
source("change/corr07/cv/cv_mse_change_sim_corr07.R")



##
## Figure 7: correlation single changepoint
##
source("change/corr07/pred/pred_mse_sim_plot.R")
source("change/corr07/coef/coef_mse_change_sim_plot.R")
source("change/corr07/cv/cv_mse_sim_plot.R")


##
## Table 5 Correlated single changepoint
##
source("change/corr07/pred/pred_mse_sim_summary.R")
source("change/corr07/coef/coef_mse_change_summary.R")
source("change/corr07/cv/cv_mse_sim_summary.R")



#########################
## Computation time
#########################
library(rbenchmark)
library("pcse")
set.seed(1999)
data("agl")
model = "within"
index = c('country', 'year')
model = 'within'
effect = 'time'

## HMBB estimation
mcmc = 1000; burn = 1000; verbose = 0; thin = 1;
formula <- growth ~ lagg1 + opengdp + openex + openimp + leftc + central + inter
##  replications = 0 causes a printing error in some machines...
a <- benchmark(
    "break 0" = {
        cp0 <- BridgeFixedPanel(formula=formula, data = agl, model = model, index = index, effect = effect,
                                mcmc=mcmc, burn=burn, verbose=verbose, Waic = TRUE, 
                                n.break = 0)},
    "break 1" = {
        cp1 <- BridgeFixedPanel(formula=formula, data = agl, model = model, index = index, effect = effect,
                                mcmc=mcmc, burn=burn, verbose=verbose, Waic = TRUE, 
                                n.break = 1)},
    "break 2" = {
        cp2 <- BridgeFixedPanel(formula=formula, data = agl, model = model, index = index, effect = effect,
                                mcmc=mcmc, burn=burn, verbose=verbose, Waic = TRUE, 
                                n.break = 2)},
    replications = 1,
    columns = c("test", "replications", "elapsed",
                "relative", "user.self", "sys.self"))

summary.table <- xtable(a, caption = "Computation Time", digits=3) 
a <- print(summary.table, booktabs = TRUE)
write(a, file = "SI/Table6.tex")

print(Sys.time() - start_time)